import numpy as np
import matplotlib.pyplot as plt
from joblib import dump, load
from sklearn.ensemble import RandomForestRegressor
import time
import pandas as pd
from sklearn.model_selection import train_test_split
import torch
device = 'cuda'
# timestamps = [[0, 20], [20, 70], [70, 169.0], [169,  200.0], [200, 230.0], [230, 291], [291, 375], [375, 500]]# 最后一个时间区间不用
# temperature0 = [175, 175, 175, 175, 175, 195, 235, 255, 255, 25, 25]
temperature1 = [173, 173, 173, 173, 173, 198, 230, 257, 257, 25, 25]
alphas = [4.29709368e-07, 6.54461698e-07, 6.03793212e-07, 5.45789093e-07, 1.01919533e-06, 6.88864878e-07, 3.23995265e-07]
hs = [7.80959263e+02, 7.37183759e+02, 7.98921888e+02, 6.94855503e+02, 8.02558178e+02, 8.30867097e+02, 7.63472858e+02]
# v = 70 / 60
H = 0.015 # milimeters
deltat = 0.5 # seconds
deltax = 0.0001 # milimeters

# 温区温度函数
def Tw(x, temperature):
    jiange = 5
    s = [25, 55.5, 60.5, 91.0, 96.0, 126.5, 131.5, 162.0, 167.0, 197.5, 202.5, 233.0, 238.0, 268.5, 273.5, 304.0, 309.0, 339.5, 344.5, 375.0, 380.0, 410.5]
    #print(s)
    for i in range(0, 22, 2):
        if x < 25:
            return x * (temperature[0] - 25) / 25
        if x >= s[i] and x <= s[i + 1]:
            return temperature[int(i/2)]
        elif x > s[i + 1] and i < 20 and x < s[i + 2]:
            return (x - s[i + 1]) / jiange * (temperature[int(i/2 + 1)] - temperature[int(i/2)]) + temperature[int(i/2)]
        
    return temperature[-1]

 
  




def model():

    df = pd.read_excel("/home/fuyusen/附件.xlsx")
    x = df['时间(s)']
    y = df['温度(ºC)']
    x_data = np.array(x)
    y_data = np.array(y)
    train, test, train_labels, test_labels = train_test_split(x_data.reshape(-1,1), y_data, test_size=0.2, random_state=42)
    rf = RandomForestRegressor(n_estimators=50, random_state=42)
    rf.fit(train, train_labels)
    #y_pred = rf.predict(test)
    score = rf.score(test, test_labels)
    # print(f"R^2 Score: {score}")
    return rf
    # import pickle
    # with open('rf.pkl', 'wb') as f:
    #     pickle.dump(rf, f, protocol=2)

def loss(alphas, hs, Td, params, flag = 0):
    temperature0 = [params[0],params[0],params[0],params[0],params[0],params[1],params[2],params[3],params[3],25,25]
    v = params[4] / 60
    timestamps = [[0, 20], [20, 70]]
    s = [25, 55.5, 60.5, 91.0, 96.0, 126.5, 131.5, 162.0, 167.0, 197.5, 202.5, 233.0, 238.0, 268.5, 273.5, 304.0, 309.0, 339.5, 344.5, 375.0, 380.0, 410.5]
    s_spot = [9, 11, 13, 17]
    for i in range(4):
        timestamps.append([timestamps[-1][1], s[s_spot[i]] / v])
    timestamps.append([timestamps[-1][1], 375])
    timestamps.append([375, 500])
    #alpha = 4.03 * 10**(-5) 
    arrt = np.arange(0, 385, deltat)
    arrx = np.arange(0, H, deltax)
    # U = np.zeros((len(arrx), len(arrt)))
    # U[:,0] = 25 # u(x_i, 0) = 25
    m = len(arrx)
    u0 = np.ones(m)
    u0 = u0 * 25
    A = np.zeros((m, m))
    A[0,0] = 0
    A[m - 1, 0] = 0
    B = np.zeros((m, m))           
    res = [25]
    res0 = [25]
    target = []
    step = 80
    B = torch.from_numpy(B).to(device)
    A = torch.from_numpy(A).to(device)
    
    # 根据alpha和h更新A, invertB
    def refresh(alpha, h):
        nonlocal A, B

        r = alpha * deltat / deltax**2
        m = B.shape[0]

        B[0,0] = 1 + h * deltax
        B[0,1] = -1
        B[m - 1, m - 1] = 1 + h * deltax
        B[m - 1, m - 2] = -1
        for i in range(1, m-1):
            A[i, i - 1] = r
            A[i, i] = 2 * (1 - r)
            A[i, i + 1] = r
        for i in range(1, m - 1):
            B[i, i - 1] = 0 - r
            B[i, i] = 2 * (1 + r)
            B[i, i + 1] = 0 - r
        B.to(device)
        A.to(device) 
        #invertB = np.linalg.inv(B)
        invertB = torch.inverse(B)
        return A, invertB
    #进行一次迭代
    def iter(u0, invertB, A, h, i):        
        def countBias(i):
            bias = np.zeros(A.shape[0])
            bias[0] = h * Tw(v * i,temperature0) * deltax
            bias[m - 1] = h * Tw(v * i,temperature0) * deltax
            bias = torch.from_numpy(bias).to(device)

            return bias

        u0 = torch.from_numpy(u0).to(device)
        bias = countBias(i)
        return invertB @ (A @ u0 + bias)
    
    in_flag = 0
    # alpha = alphas[0]
    # h = hs[0]
    # A, invertB = refresh(alpha, h)
    for i , curr_time in enumerate(arrt):
        if i == len(arrt) - 1:
            break
        if curr_time >= timestamps[flag][1]:
            break
        if curr_time >= timestamps[in_flag][0]:
            # print(f"refresh:{in_flag + 1}")
            if in_flag > flag:
                break
            alpha = alphas[in_flag]
            h = hs[in_flag]
            A, invertB = refresh(alpha, h)
            in_flag = in_flag + 1
            #print(arrt[i])
        iter_result = iter(u0, invertB, A, h, arrt[i + 1])
        iter_result = iter_result.cpu().numpy()
        
        # if int(len(arrx) / 2) > boundaries[flag]:
        #     refresh(alpha[flag], alpha[flag])
        #print(result)
        
        res.append(u0[int(len(arrx) / 2)])
        res0.append(u0[0])
        u0 = iter_result
        
        if arrt[i] >= 19 :
            target.append(Td(arrt[i]))
            # if in_flag == flag + 1:
            #     if i % step == 0:
            #         # print("exp: " + str(res[-1]))
            #         # print("exp_0: " + str(res0[-1]))
            #     if i % step == 0 :
            #         # print("standard:" + str(target[-1]))

    res = res[-len(target):]    
    return res

rf = model()
def Td(t):
    return rf.predict([[t]])[0]

# %%
# 求取面积的函数
def S(temps):
    sum = 0
    for i in range(len(temps)):
        sum += (temps[i] - 217) * 0.5
    
    return sum

def Symmetryloss(temps,tstart,tmax,tend):
    sum1 = 0
    for i in range(tstart+1,tmax):
        sum1 += (temps[i]-217) * 0.5

    S1 = sum1 + ((temps[tstart]+temps[tmax]) / 2 - 217) * 0.5

    sum2 = 0
    for i in range(tmax+1,tend):
        sum2 += (temps[i] - 217) * 0.5

    S2 = sum2 + ((temps[tmax]+temps[tend]) / 2 - 217) * 0.5

    operand1 = abs(S1-S2) / max(S1,S2)
    operand2 = abs(tend + tstart - 2 * tmax) / max(tmax-tstart,tend-tmax)

    sigma = max(operand1,operand2)

    return sigma
    
def getSection(temps):
    temp_max = max(temps)
    start = 0
    end = 0
    for i in range(len(temps)):
        if temps[i] >= 217 and start == 0:
            start = i

        if temps[i] == temp_max:
            end = i
            break

    return temps[start:end+1]

# %%
import pandas as pd
import numpy as np
def check(time,temps):
    time_150 = 0
    # time_150_found = False
    time_190 = 0
    # time_190_found = False
    time_217_start = 0
    tstart = 0
    # time_217_start_found = False
    time_217_end = 0
    tend = 0
    # time_217_end_found = False
    max_temp = max(temps)
    tmax = temps.index(max_temp)
    for i in range(len(time)):
        if (temps[i] >= 150 and time_150 == 0):
            time_150 = time[i]
            

        if (temps[i] >= 190 and time_190 == 0):
            time_190 = time[i]

        if (temps[i] >= 217 and time_217_start == 0):
            time_217_start = time[i]
            tstart = i

        if (temps[i] < 217 and time_217_start != 0):
            time_217_end = time[i]
            tend = i-1
            break

    time_150to190 = time_190 - time_150
    time_over217 = time_217_end - time_217_start
    grad = np.gradient(temps) * 2
    grad_check = (abs(grad) <= 3).all()
    print("gradient",grad_check)
    time_check1 = (time_150to190 >= 60 and time_150to190 <= 120)
    print("150-190",time_150to190,time_check1)
    time_check2 = (time_over217 >= 40 and time_over217 <= 90)
    print("over217",time_over217,time_check2)
    temp_check = (max_temp >= 240 and max_temp <= 250)
    print("maxtemp",max_temp,temp_check)
    symmetry = Symmetryloss(temps,tstart,tmax,tend)
    symmetry_check = (symmetry <= 0.3714286*1.05)
    print("对称性：",symmetry,symmetry_check)
    res = grad_check and time_check1 and time_check2 and temp_check and symmetry_check
    # print("检查完毕！")
    return res



        

# %%
import math
import random 
import matplotlib.pyplot as plt

# def func(x, y):                  #函数优化问题
#     res= 4*x**2-2.1*x**4+x**6/3+x*y-4*y**2+4*y**4
#     return res
#x为公式里的x1,y为公式里面的x2
class SA:
    def __init__(self, func=S, getTemps=loss, iter=10, T0=100, Tf=0.1, alpha=0.94):
        self.func = S
        self.getTemps = loss
        self.iter = iter         #内循环迭代次数,即为L =100
        self.alpha = alpha       #降温系数，alpha=0.99
        self.T0 = T0             #初始温度T0为100
        self.Tf = Tf             #温度终值Tf为0.01
        self.T = T0              #当前温度
        # self.x = [random() * 11 -5  for i in range(iter)] #随机生成100个x的值
        # self.y = [random() * 11 -5  for i in range(iter)] #随机生成100个y的值
        self.temp_zone1_5 = [ random.uniform(165,185) for i in range(iter)]
        self.temp_zone6 = [random.uniform(185,205) for i in range(iter)]
        self.temp_zone7 = [random.uniform(225,245) for i in range(iter)]
        self.temp_zone8_9 = [random.uniform(245,265) for i in range(iter)]
        self.speed = [random.uniform(65,100) for i in range(iter)]
        self.most_best =[]
        """
        random()这个函数取0到1之间的小数
        如果你要取0-10之间的整数（包括0和10）就写成 (int)random()*11就可以了，11乘以零点多的数最大是10点多，最小是0点多
        该实例中x1和x2的绝对值不超过5（包含整数5和-5），（random() * 11 -5）的结果是-6到6之间的任意值（不包括-6和6）
        （random() * 10 -5）的结果是-5到5之间的任意值（不包括-5和5），所有先乘以11，取-6到6之间的值，产生新解过程中，用一个if条件语句把-5到5之间（包括整数5和-5）的筛选出来。
        """
        self.history = {'f': [], 'T': [], 'idx':[]}

    def generate_new(self, temp_zone1_5,temp_zone6, temp_zone7, temp_zone8_9,speed):   #扰动产生新解的过程
        count = 0
        while count <= 5:
            # x_new = x + self.T * (random() - random())
            # y_new = y + self.T * (random() - random())
            temp_zone1_5_new = temp_zone1_5 + 3 * random.random() - 1
            temp_zone6_new = temp_zone6 + 3 * random.random() - 1
            temp_zone7_new = temp_zone7 + 3 * random.random() - 1
            temp_zone8_9_new = temp_zone8_9 + 3 * random.random() - 1
            speed_new = speed + 3 * random.random() - 1
            # print("扰动完毕！")
            if (165 <= temp_zone1_5_new <= 185) and (185 <= temp_zone6_new <= 205) and (225 <= temp_zone7_new <= 245) and (245 <= temp_zone8_9_new <= 265) and (65 <= speed <= 100):  
                # print(temp_zone1_5_new, temp_zone6_new, temp_zone7_new, temp_zone8_9_new, speed_new)
                temps = self.getTemps(alphas,hs, Td,[temp_zone1_5_new, temp_zone6_new, temp_zone7_new, temp_zone8_9_new, speed_new],6)
                # print("获取用于检查制程界限的温度完毕！")
                time = np.linspace(0,375,len(temps))
                # plt.plot(time,temps)
                # plt.show()
                # print(len(temps),len(time))
                if check(time,temps):
                    return temp_zone1_5_new, temp_zone6_new, temp_zone7_new, temp_zone8_9_new, speed_new     
                print("不满足制程条件！")
            count += 1                             #重复得到新解，直到产生的新解满足约束条件
        return 176.75,186.86, 227.58, 264.52, 88

    def Metrospolis(self, f, f_new):   #Metropolis准则
        # if not check(temps):
        #     return 0
        print("旧面积：",f, "新面积：",f_new)
        if f_new <= f:
            print("M检查完毕，采用新解")
            return True
        else:
            p = math.exp((f - f_new) / self.T)
            if random.random() < p:
                print("M检查完毕，采用新解")
                return True
            else:
                print("M检查完毕，不采用新解")
                return False

    def best(self,f_list):    #获取最优目标函数值
        # f_list = []    #f_list数组保存每次迭代之后的值
        # for i in range(self.iter):
        #     temps = self.getTemps(alphas,hs,Td,[self.temp_zone1_5[i],self.temp_zone6[i],self.temp_zone7[i],self.temp_zone8_9[i],self.speed[i]],7)
        #     temps_section = getSection(temps)

        #     f = self.func(temps_section)
        #     f_list.append(f)
        f_best = min(f_list)
        
        idx = f_list.index(f_best)
        return f_best, idx    #f_best,idx分别为在该温度下，迭代L次之后目标函数的最优解和最优解的下标

    def run(self):
        count = 0
        #外循环迭代，当前温度小于终止温度的阈值
        while self.T > self.Tf:       
            print("外循环第",count,"轮")
            f_list = []
            #内循环迭代100次
            for i in range(self.iter): 
                print("内循环第",i,"轮")
                temps_ori = self.getTemps(alphas,hs,Td,[self.temp_zone1_5[i],self.temp_zone6[i],self.temp_zone7[i],self.temp_zone8_9[i],self.speed[i]],6)
                f = 0
                if count == 0:
                    if not check(np.linspace(0,375,712),temps_ori):
                        f = 550

                    else:
                        temps_ori_sec = getSection(temps_ori)
                        f = self.func(temps_ori_sec)

                else:
                    temps_ori_sec = getSection(temps_ori)
                    f = self.func(temps_ori_sec)
                # print("旧温度切割完毕！")
                temp_zone1_5_new, temp_zone6_new, temp_zone7_new, temp_zone8_9_new, speed_new = self.generate_new(self.temp_zone1_5[i],self.temp_zone6[i],self.temp_zone7[i],self.temp_zone8_9[i],self.speed[i]) #产生新解
                temps = self.getTemps(alphas,hs,Td,[temp_zone1_5_new, temp_zone6_new, temp_zone7_new, temp_zone8_9_new, speed_new],6)
                # print("新温度获取完毕！")
                temps_sec = getSection(temps)
                f_new = self.func(temps_sec)                        #产生新值
                if self.Metrospolis(f, f_new):                         #判断是否接受新值
                    self.temp_zone1_5[i] = temp_zone1_5_new           #如果接受新值，则把新值的x,y存入x数组和y数组
                    self.temp_zone6[i] = temp_zone6_new
                    self.temp_zone7[i] = temp_zone7_new
                    self.temp_zone8_9[i] = temp_zone8_9_new
                    self.speed[i] = speed_new
                    f_list.append(f_new)

                else:
                    f_list.append(f)

            # 迭代L次记录在该温度下最优解
            ft, idx = self.best(f_list)
            print("该温度下面积最小值：",ft)
            print("该温度下最优解：",sa.temp_zone1_5[idx],sa.temp_zone6[idx],sa.temp_zone7[idx],sa.temp_zone8_9[idx],sa.speed[idx])
            self.history['f'].append(ft)
            self.history['T'].append(self.T)
            self.history['idx'].append(idx)
            #温度按照一定的比例下降（冷却）
            self.T = self.T * self.alpha        
            count += 1
        

            # 得到最优解
        # f_best, idx = self.best()
        # print(f"F={f_best}, temp_zone1_5={self.temp_zone1_5[idx]}, temp_zone6={self.temp_zone6[idx]},temp_zone7={self.temp_zone7[idx]}, temp_zone8_9={self.temp_zone8_9[idx]}, speed={self.speed[idx]}")

        


# %%
sa = SA()
sa.run()
idx = sa.history["idx"][-1]
print("最小面积：",sa.history["f"][-1])
print("对应变量值：")
print("温区1-5温度：",sa.temp_zone1_5[idx])
print("温区6温度：",sa.temp_zone6[idx])
print("温区7温度：",sa.temp_zone7[idx])
print("温区8-9温度：",sa.temp_zone8_9[idx])
print("过炉速度：",sa.speed[idx])